/*
 
 Trabalho 1
 Programação Concorrente
 
 Cálculo de PI com precisão de 10 milhões de dígitos
 Método: Borwein (Sequencial)
 Biblioteca de precisão: The GNU Multiple Precision Arithmetic Library (GMP)
 
 */


#include <stdio.h>
#include <gmp.h>

#define LOG_2_10 3.32192809488736234787

int main (void)
{
    /* Contador de iterações */
    int i = 1;
    /* Fator multiplicativo */
    int k = 8;
    /* Vatiáveis big float da GMP */
    mpf_t pi_old, pi_new;
    mpf_t a_old, a_new;
    mpf_t y_old, y_new;
    mpf_t temp1, temp2;
    mpf_t epsilon;
    /* Variável para gravar o nome do arquivo de saída */
    char filename[10];
    /* Ponteiro para arquivo de saída */
    FILE * F;
    
    /*
     
     Define a precisão padrão usada nos cálculos:
     
     OBS: Utilizamos 10000100 casas de precisão para
     possíveis overflows durante o cálculo.
     Como o argumento é em bits, e não em casas decimais,
     multiplicamos por uma constante log2(10),
     que é a quantidade média de bits por casa decimal
     
     */
    mpf_set_default_prec(LOG_2_10 * 10000100);
    
    /*
     
     Inicializa as variáveis de acordo com o método:
     a = 6-4*sqrt(2)
     y = sqrt(2)-1
     
     */
    
    mpf_init(pi_old);
    mpf_init(pi_new);
    mpf_init(a_old);
    mpf_init(y_old);
    mpf_init(temp1);
    mpf_init(temp2);
    
    mpf_init_set_d(a_new, 32.0);
    mpf_sqrt(a_new, a_new);
    mpf_ui_sub(a_new, 6, a_new);
    
    mpf_init_set_d(y_new, 2.0);
    mpf_sqrt(y_new, y_new);
    mpf_sub_ui(y_new, y_new, 1);
    
    mpf_init_set_str(epsilon, "1e-10000000", 10);
    
    mpf_ui_div(pi_new, 1, a_new);
    
    /*
     
     Calcula as demais iterações, gravando cada resultado parcial
     em um arquivo de nome it-x.txt
     
     */
    do
    {
        mpf_swap(pi_old, pi_new);
        mpf_swap(a_old, a_new);
        mpf_swap(y_old, y_new);
       
        mpf_pow_ui(y_old, y_old, 4);
        mpf_ui_sub(y_old, 1, y_old);
        mpf_sqrt(y_old, y_old);
        mpf_sqrt(y_old, y_old);
        mpf_add_ui(temp1, y_old, 1);
        mpf_ui_sub(y_new, 1, y_old);
        mpf_div(y_new, y_new, temp1);
        
        mpf_add_ui(temp1, y_new, 1);
        
        mpf_pow_ui(temp2, y_new, 2);
        mpf_add(temp2, temp2, temp1);
        mpf_mul(temp2, temp2, y_new);
        mpf_mul_ui(temp2, temp2, k);
        k *= 4;
        
        mpf_pow_ui(temp1, temp1, 4);
        mpf_mul(temp1, temp1, a_old);
        mpf_sub(a_new, temp1, temp2);
        
        mpf_ui_div(pi_new, 1, a_new);
        
        mpf_sub(pi_old, pi_new, pi_old);
        mpf_abs(pi_old, pi_old);
        sprintf(filename, "it-%d.txt", i++);
        F = fopen(filename, "w");
        gmp_fprintf(F, "%.10000000Ff\n", pi_new);
        fclose(F);
    } while ( mpf_cmp(epsilon, pi_old) < 0 );
    
    /* Libera o espaço de memória alocado para as estruturas */
    mpf_clear(pi_old);
    mpf_clear(pi_new);
    mpf_clear(a_old);
    mpf_clear(a_new);
    mpf_clear(y_old);
    mpf_clear(y_new);
    mpf_clear(temp1);
    mpf_clear(temp2);
    mpf_clear(epsilon);
    
    return 0;
}
